home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / exp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  8.5 KB  |  364 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  * A replacement exp() routine for PML library.  An original algorithm  *
  4.  * is not changed but a table of fractionals power of 2 was recalcul-   *
  5.  * ated (believe it or not - this has an influence on a final result).  *
  6.  * Also divisions by power of 2 were replaced by applications of ldexp  *
  7.  * routine which is faster and better preserves precision               *
  8.  *                                                                      *
  9.  * Michal Jaegermann, December 1990                                     *
  10.  *                                                                      *
  11.  ************************************************************************
  12.  *                                                                      *
  13.  * 68881 support added by Michael Ritzert                *
  14.  *                                                                      *
  15.  ************************************************************************
  16.  */
  17.   /************************************************************************
  18.  *                                    *
  19.  *                N O T I C E                *
  20.  *                                    *
  21.  *            Copyright Abandoned, 1987, Fred Fish        *
  22.  *                                    *
  23.  *    This previously copyrighted work has been placed into the    *
  24.  *    public domain by the author (Fred Fish) and may be freely used    *
  25.  *    for any purpose, private or commercial.  I would appreciate    *
  26.  *    it, as a courtesy, if this notice is left in all copies and    *
  27.  *    derivative works.  Thank you, and enjoy...            *
  28.  *                                    *
  29.  *    The author makes no warranty of any kind with respect to this    *
  30.  *    product and explicitly disclaims any implied warranties of    *
  31.  *    merchantability or fitness for any particular purpose.        *
  32.  *                                    *
  33.  ************************************************************************
  34.  */
  35.  
  36.  
  37. /*
  38.  *  FUNCTION
  39.  *
  40.  *    exp   double precision exponential
  41.  *
  42.  *  KEY WORDS
  43.  *
  44.  *    exp
  45.  *    machine independent routines
  46.  *    math libraries
  47.  *
  48.  *  DESCRIPTION
  49.  *
  50.  *    Returns double precision exponential of double precision
  51.  *    floating point number.
  52.  *
  53.  *  USAGE
  54.  *
  55.  *    double exp (x)
  56.  *    double x;
  57.  *
  58.  *  REFERENCES
  59.  *
  60.  *    Fortran IV plus user's guide, Digital Equipment Corp. pp B-3
  61.  *
  62.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  63.  *    1968, pp. 96-104.
  64.  *
  65.  *  RESTRICTIONS
  66.  *
  67.  *    Inputs greater than log(MAXDOUBLE) result in overflow.
  68.  *    Inputs less than log(MINDOUBLE) result in underflow.
  69.  *
  70.  *    The maximum relative error for the approximating polynomial
  71.  *    is 10**(-16.4).  However, this assumes exact arithmetic
  72.  *    in the polynomial evaluation.  Additional rounding and
  73.  *    truncation errors may occur as the argument is reduced
  74.  *    to the range over which the polynomial approximation
  75.  *    is valid, and as the polynomial is evaluated using
  76.  *    finite precision arithmetic.
  77.  *    
  78.  *  PROGRAMMER
  79.  *
  80.  *    Fred Fish
  81.  *
  82.  *  INTERNALS
  83.  *
  84.  *    Computes exponential from:
  85.  *
  86.  *        exp(x)    =    2**y  *  2**z  *  2**w
  87.  *
  88.  *    Where:
  89.  *
  90.  *        y    =    int ( x * log2(e) )
  91.  *
  92.  *        v    =    16 * frac ( x * log2(e))
  93.  *
  94.  *        z    =    (1/16) * int (v)
  95.  *
  96.  *        w    =    (1/16) * frac (v)
  97.  *
  98.  *    Note that:
  99.  *
  100.  *        0 =< v < 16
  101.  *
  102.  *        z = {0, 1/16, 2/16, ...15/16}
  103.  *
  104.  *        0 =< w < 1/16
  105.  *
  106.  *    Then:
  107.  *
  108.  *        2**z is looked up in a table of 2**0, 2**1/16, ...
  109.  *
  110.  *        2**w is computed from an approximation:
  111.  *
  112.  *            2**w    =  (A + B) / (A - B)
  113.  *
  114.  *            A    =  C + (D * w * w)
  115.  *
  116.  *            B    =  w * (E + (F * w * w))
  117.  *
  118.  *            C    =  20.8137711965230361973
  119.  *
  120.  *            D    =  1.0
  121.  *
  122.  *            E    =  7.2135034108448192083
  123.  *
  124.  *            F    =  0.057761135831801928
  125.  *
  126.  *        Coefficients are from HART, table #1121, pg 206.
  127.  *
  128.  *        Effective multiplication by 2**y is done by a
  129.  *        floating point scale with y as scale argument.
  130.  *
  131.  */
  132.  
  133. #if !defined (__M68881__) && !defined (sfp004)    /* mjr++        */
  134.  
  135. #include <stdio.h>
  136. #include <math.h>
  137. #include "pml.h"
  138.  
  139. static char funcname[] = "exp";
  140.  
  141. # define C  20.8137711965230361973    /* Polynomial approx coeff.    */
  142. /* # define D  1.0    */              /* Polynomial approx coeff.    */
  143. # define E  7.2135034108448192083    /* Polynomial approx coeff.    */
  144. # define F  0.057761135831801928    /* Polynomial approx coeff.    */
  145.  
  146.  
  147. /************************************************************************
  148.  *                                    *
  149.  *    This table is fixed in size and reasonably hardware        *
  150.  *    independent.  The given constants were computed on Atari ST     *
  151.  *      using integer arithmetic and decimal values for 2**(1/2),       *
  152.  *      2**(1/4) and 2**(1/8) taken from Appendix C of HART et al.      *
  153.  *                                    *
  154.  ************************************************************************
  155.  */
  156.  
  157. static double fpof2tbl[] = {
  158.     1.0/*0000000000000000000*/,        /*    2 ** 0/16    */
  159.     1.04427378242741384032,        /*    2 ** 1/16    */
  160.     1.09050773266525765920,        /*    2 ** 2/16    */
  161.     1.13878863475669165369,        /*    2 ** 3/16    */
  162.     1.18920711500272106671,        /*    2 ** 4/16    */
  163.     1.24185781207348404858,        /*    2 ** 5/16    */
  164.     1.29683955465100966592,        /*    2 ** 6/16    */
  165.     1.35425554693689272828,        /*    2 ** 7/16    */
  166.     1.41421356237309504880,        /*    2 ** 8/16    */
  167.     1.47682614593949931138,        /*    2 ** 9/16    */
  168.     1.54221082540794082360,        /*    2 ** 10/16    */
  169.     1.61049033194925430816,        /*    2 ** 11/16    */
  170.     1.68179283050742908605,        /*    2 ** 12/16    */
  171.     1.75625216037329948310,        /*    2 ** 13/16    */
  172.     1.83400808640934246346,        /*    2 ** 14/16    */
  173.     1.91520656139714729384        /*    2 ** 15/16    */
  174. };
  175.  
  176.     char * xcptstr;
  177.     struct exception xcpt;
  178.  
  179. double exp (x)
  180. double x;
  181. {
  182.     register int index;
  183.     int y;
  184.     double w;
  185.     double a;
  186.     double b;
  187.     double temp;
  188.  
  189.     if (x > LOGE_MAXDOUBLE) {
  190.     xcpt.type = OVERFLOW;
  191.     xcpt.retval = MAXDOUBLE;
  192.     xcptstr = "OVER";
  193.     goto eset;
  194.     }
  195.     if (x <= LOGE_MINDOUBLE) {
  196.     xcpt.type = UNDERFLOW;
  197.     xcpt.retval = 0.0;
  198.     xcptstr = "UNDER";
  199. eset:
  200.     xcpt.name = funcname;
  201.     xcpt.arg1 = x;
  202.     if (!matherr (&xcpt)) {
  203.             errno = ERANGE;
  204.         fprintf (stderr, "%s: %sFLOW error\n", funcname, xcptstr);
  205.     }
  206.     } else {
  207.     x *= LOG2E;
  208.     w = ldexp(modf(x,&a), 4);
  209.     y = a;
  210.     w = ldexp(modf(w, &temp), -4);
  211.     index = temp;
  212.     
  213.     b = w * w;
  214.     a = C + b;
  215.     b = w * (E + (F * b));
  216.     xcpt.retval = ldexp (((a + b) / (a - b)) *
  217.         (index < 0 ? ldexp(fpof2tbl[16 + index], -1) : fpof2tbl[index]), y);
  218.     }
  219.     return (xcpt.retval);
  220. }
  221. #endif    /* !__M68881__ && !sfp004    */
  222. #ifdef    sfp004
  223.  
  224. __asm("
  225.  
  226. comm =     -6
  227. resp =    -16
  228. zahl =      0
  229.  
  230. ");    /* end asm    */
  231.  
  232. #endif    sfp004
  233. #if defined (__M68881__) || defined (sfp004)
  234.     __asm(".text; .even");
  235.  
  236. # ifdef    ERROR_CHECK
  237.  
  238.     __asm("
  239.  
  240.  
  241. _Overflow:
  242.     .ascii \"OVERFLOW\\0\"
  243. _Domain:
  244.     .ascii \"DOMAIN\\0\"
  245. _Error_String:
  246.     .ascii \"exp: %s error\\n\\0\"
  247. .even
  248.  
  249. | pml compatible expgent
  250. | m.ritzert 7.12.1991
  251. | ritzert@dfg.dbp.de
  252. |
  253. |    /* NAN  = {7fffffff,ffffffff}        */
  254. |    /* +Inf = {7ff00000,00000000}        */
  255. |    /* -Inf = {fff00000,00000000}        */
  256. |    /* MAX_D= {7fee42d1,30773b76}        */
  257. |    /* MIN_D= {ffee42d1,30773b76}        */
  258.  
  259. .even
  260. double_max:
  261.     .long    0x7fee42d1
  262.     .long    0x30273b76
  263. double_min:
  264.     .long    0xffee42d1
  265.     .long    0x30273b76
  266. NaN:
  267.     .long    0x7fffffff
  268.     .long    0xffffffff
  269. p_Inf:
  270.     .long    0x7ff00000
  271.     .long    0x00000000
  272. m_Inf:
  273.     .long    0xfff00000
  274.     .long    0x00000000
  275. ");
  276. #endif    ERROR_CHECK
  277.  
  278. __asm("
  279. .even
  280.     .globl _exp
  281. _exp:
  282. ");    /* end asm    */
  283.  
  284. #endif    /* __M68881__ || sfp004    */
  285. #ifdef    __M68881__
  286.  
  287.     __asm("
  288.     fetoxd    a7@(4), fp0    | exp
  289.     fmoved    fp0,a7@-    | push result
  290.     moveml    a7@+,d0-d1    | return_value
  291.     ");    /* end asm    */
  292.  
  293. #endif    __M68881__
  294. #ifdef    sfp004
  295.     __asm("
  296.     lea    0xfffa50,a0
  297.     movew    #0x5410,a0@(comm)    | specify function
  298.     cmpiw    #0x8900,a0@(resp)    | check
  299.     movel    a7@(4),a0@        | load arg_hi
  300.     movel    a7@(8),a0@        | load arg_low
  301.     movew    #0x7400,a0@(comm)    | result to d0
  302.     .long    0x0c688900, 0xfff067f8    | wait
  303.     movel    a0@,d0
  304.     movel    a0@,d1
  305.     ");    /* end asm    */
  306.  
  307. #endif    sfp004
  308. #if defined (__M68881__) || defined (sfp004)
  309. # ifdef    ERROR_CHECK
  310.     __asm("
  311.     lea    double_max,a0    |
  312.     swap    d0        | exponent into lower word
  313.     cmpw    a0@(16),d0    | == NaN ?
  314.     beq    error_nan    |
  315.     cmpw    a0@(24),d0    | == + Infinity ?
  316.     beq    error_plus    |
  317.     swap    d0        | result ok,
  318.     rts            | restore d0
  319. ");
  320. #ifndef    __MSHORT__
  321. __asm("
  322. error_plus:
  323.     swap    d0
  324.     moveml    d0-d1,a7@-
  325.     movel    #63,_errno    | NAN => errno = EDOM
  326.     pea    _Overflow    | for printf
  327.     bra    error_exit    |
  328. error_nan:
  329.     moveml    a0@(24),d0-d1    | result = +inf
  330.     moveml    d0-d1,a7@-
  331.     movel    #62,_errno    | NAN => errno = EDOM
  332.     pea    _Domain        | for printf
  333. ");
  334. #else    __MSHORT__
  335. __asm("
  336. error_plus:
  337.     swap    d0
  338.     moveml    d0-d1,a7@-
  339.     movew    #63,_errno    | NAN => errno = EDOM
  340.     pea    _Overflow    | for printf
  341.     bra    error_exit    |
  342. error_nan:
  343.     moveml    a0@(24),d0-d1    | result = +inf
  344.     moveml    d0-d1,a7@-
  345.     movew    #62,_errno    | NAN => errno = EDOM
  346.     pea    _Domain        | for printf
  347. ");
  348. #endif    __MSHORT__
  349. __asm("
  350. error_exit:
  351.     pea    _Error_String    |
  352.     pea    __iob+52    |
  353.     jbsr    _fprintf    |
  354.     addl    #12,a7        |
  355.     moveml    a7@+,d0-d1
  356.     rts
  357.     ");
  358. # else    ERROR_CHECK
  359.  
  360. __asm("rts");
  361.  
  362. # endif    ERROR_CHECK
  363. #endif /* __M68881__ || sfp004    */
  364.